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IN TRODUCTION 

This report contains the documentation for the LPASS program. 
It conalsta,of- the design procedure used, a description of the 
program, and deaign examples using the program. • 

.The purpose of the LPASS program is the design of. a maximally 
flat flutterwtorth or an equiripple Chebychev lowpass d*lgital filter. 
Starting with an analog filter, the bilinear Z t.ransform is used to ' 
f ind '^n^equivalent digital filter. The user enters the following 
parameter^: the nimber of second order sections, the type of filter, 

the sampling interval, the -3db cutoff frequency, the*starting frequency 

' ' ■ . ■ . ■ ■> 

and the frequency incremient. H a Chebychev filter is being designed, 
the ripple must also be entered. 

The program Cj|lculates the digital filter coefficients for up to 
three second order sections in cascade. The program is designed to 
calculate up to a sixth order filter, thus. the filter order' is two 

* ■ • 

tiiaes the number of cascaded second oirA^r sections. The -filter 
magnitude response is generated over thm frequency ^nteryal specified 
by the input* 

The LPAS^ program, written in Fortran IV, is supplied as a card 
deck with this 1^ report. ^ The program is in the form of a subroutine 
and can be used as is by a call statement from the main prog'ram. 
Data i;iay be input via cards with output available through a lint?^ 
printerf The input/output devices may be altered as^explained in 
this report. Graphics routines may easily be appended to the program. 



I. Dfesign Procedure 

J — 

* A. Prelfmlnary Discussion 

The transfer function of a second order digital filter In the 

Z domain Is given by ^ 

K,(ArtZ'^+A,Z+A,) 

H(Z) - (1) 

Z +Bj^Z+B2 

where the A's and B's are the coefficients of the numerator and 
denominator respectively. One conanon method of designing a digital 
filter is to start with an analog transfer function H(S) and 
transform it to the digital transfer function 'H{Z>. This 
program will calculate the scale factor and the coefficients 

Aq, A^, A2, Bj^, and B^. The transformation used is the extended 
bilinear Z transform defined as 

^ T ^Z+1^ * • 

.where T, is t.ht .sampling internal. VHien this transform is employed, 
the desired frequencies must first be prewarjped to iooake them com- 
patible with the digital filter. The prewarped cutoff frequency is 

. given by • . ■ 

0) T 

mC ^ - Tan (-r|-) •. ' (3) 

This prewarping is clone by the program. 
B. Butterworth Low-Pass FIX tar 

We start with a normalized second order low-pass filter in 
the S plane. ' f 



H(S) - -2 ■ 

S + 2Scose +1 • ' 



6 



where the angle 9 Is in degrees (in the pro^gram) . 6 may be found 



from the Butterworth circle and the relationship 

_ ' +jip<2m-l)/2n 



(5) ^ 



where n is the order of the filter and m « !» 2, 3^ • . • , n. 

* ' * »^ * 

This relationship is determined by the following procedure. By 

definition, a filter is n order Butterworth loift-pass ,if its gain 
characteristic is ^ . * " 



(6) 



1> (-^) 
c 



2n 



where a Is the gain, o) is the desired cutoff frequency i^and n is 

. ■ . ' c - ■ ' • • • 

the order of the filter. Note that |H^(ja>)| goes to zero as o) goes 
to Infinity, indicating the filter does attenuate the higher frequencies. 
To determine its efficiency as a low-pas^ f ilte^ we calculate 



an 



c 



2n-l 



l+(~) 

OS 

c 



In! 3/2 



(7) 



Thus 



_d. 



H„(ju))1 



n 



0 



(8) 



for all n and hence the gain characteristic stays flat for u) close 
to 0. Also 



7^|h (jo)) 



an 



20) /2 
c 



(9) 



and hence, the decline rate Or "roll-off "-of the gain characteristic 



at w - c# becomes sharper as n increases. In other words, the 

■c 

approximation to the ideal Xow^pass filter isproves for larger n. 
The order n is qhosen according to desired specifications. The ^ 
references have equations, curves, and tables that select n. given 
the spec if ligations. Far example, page 227 of Rabiner and Gold gives 

i 

an equation for calculating n when the transition band is specified. 

In the design, the poles for the full frequency response, H(S), 
of the n*^*^ order Butterworth filter mist be determined. The pro- 
cedure is as follows: . a 



[H(S)H(-S)] 



2n 

1 + (-^) J 



a 



2n 



S 1 + (— ) 



a 













n. 


1 + - 








2 






Hi 




c 





r 2-1 
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n 


1 + 








- -2- 






CO 






1- c . 




f 



a 



1 - 



^ c 



, for n even 



for. n odd 



(10) 



Setting the denominators equal to zero, 



0) 



(11) 



Thus, the pole locations are the 2n roots of ±1, depending on whether 
the order is odd or even. These roots are located op a citcle with ■ 
radius w centered at the otigin of the S plane arid have symmetry 



ERIC 



6^ 



with respect both real and Imaginary axes. For n odd, a pair of ^ 
- roots are on the real axis and the ycat are separated by ir/n radians. 
Fot n even* a pair of roots are located tr/2n radians from t^e real 
axis and the rest are again separated by ir/n radians. No roots are 
on the imaginary axis for either even or odd n. 

Let p, , . . . , Pj he the roots. From the syiametry of the pole 

locations, if p^ p^^ are the roots lying in the right-half 

plane, the left-half plane roots are -p^^, .... -p^r Th® 
oagnitude-squared function can then he written as 

a <-l) 0)^2. 

H^(S)H^(-S) - (s^^j ^ ^ ^ (S+p^XS-p^) . . . (S^) • 
To be stable, H^(S) must have all its poles In, the left-half plane, thu 



n 
c 



aw_ . ' (13) 



"n^^^ " (S+p^) . . . (S+p^) > 
The program is written with unify gain at D6^(<ij«0), • therefore a = 1. 

In order to locate the poles as specified^ above, consider the 
following set of equations. ' 

I „ _g±jTr(2in-l) m « 1, 2, . . . , n; for n even 

(U) 

.1 » „e*^ , k - 0, 1, . . . , n; for n odd 

Substituting equations (14) into equation (11) yields 

' " - r S, ±j7T(2m-l)/2n . s ' ^ * 

f — 1' « -e ^ , tn * 1, 2, . . . , n; for n even 

. (15) 

r^l - ^g±jiTk/n ic = 0, 1, . , „ , n; for n odd- ' 
c • •■ 

Equations (15) will give the pole locations des$:ribed, d^ove^ 

Consider the fona of equations (15) ^ 

S - -oj e*^® - 0) [-COS0 i'jsinO]. - . (16) 

c • c , 



From this relationship, it cao be seen that the magnitude ^or each pole 
is u) , regardless of the angle, and thus all the* poles lie cm a circle 



c 

with radius o) « 
c 



As an example consider a second o^der filter, n ■ 2, ^ 

P ^±5Tr(2m-l)/4 m « 1, 2 

c 



c 



tl c — ^ 



S+2 " /^135 ' 



The relationship of these roots about the circle of radium w 
is ^lustrated id Figure 1. The angle 8 is always joeasured from 
the negative real axis. , ^ , 

In the program, only the angle(s) less than 90** are considered * 
so that poles lie in the left-hal£ plane since poles in the left-half 

r ' ' ' -K 

ft ' ' * 

plane are stable. Putting .8 - 45° into equation (4) yields poles 
at -0.70.7 ± j0.707. These locations are in the ieft-half plane. 

In the program, only even drder filters are considered. 

Below are the values of 6 for 1, 2, and 3 second. order sections 
In cascade. " 



Cascaded Filter 
Sections Order 



Angle 

N . " n eV 



.1 2 45'. 

2 4^^ .22.5% 67.5** 

a 6 75% 45% 15 



O 



These calculated angles a^e incorporated in the program in tf^e order 

■ • - . 

given above. * >. 



For N second drder sections there are N 9's. Only one specific 

i ' • , ■ 

Q is us(?d per stage/ l^cause each stage has only one set of pole ^ 

lacatioi^s.t ' 

* ■ • ' ' ■ . ^ ' "Ml 

Thfe' following is the procedure to derive the saagnitude of the 1 - 

■ ' /. ' ' ' ■ ' 

stage, where i varies froia 1 to N.. 

Given the nonoal^zed second order low-pass transfer function 

equation (A), we employ" the low-pass to low-pass transfonoation for 

an arbitrary cutoff frequency o) given by . . * ; 

. S -® • * (17 ) 

c 

For. the i &tage,% equation (4) becomes * 

2 . • ^ ' 

KAS) ' 2 - (18) 

, ' , S +2SU1 cose .-hu • . 

c 1 c . 

The extended bilinear Z transform, equation (2), Is used to get to 

the digital domain. Employing equation (2)^ on equation (18) and sub- 

stitutlngiroc for w yields ' 

^' -~(Z^-2Z-I-X)-4(Z -1)WDC cose +WDC- (2 +2Z+1) 

Putting the denominator ^ equation (19) in monlc form yields the 
transfer function for the i^ stage of the filter 

H. (Z) - -iV^ . . .(20) 

Z+B^,Z+B2, 

Equation (20) is the san» as equation (1) with the exception of the 
subscripts. Inequation (20) 



- -4 + 4 WDC cos8 . + 'WDC^ (21) 
1 ^2 Tsv i . 

2 



»2i 





8 

■ 






4 4 


WDCcosS^+WDC^ 





ST 

Letting Z - a and S « jw and taking the magnitude of H^(Jw) we 
have 

/(AqCOS (2ajT)+Aj^cos (wtHA2 ) ^+ (AqS in (2uT)+A^8int^«»T) ) ^ 
J^^'cos (wD+Bg^) ^+(Sin (2wT)+Bj,^^ 



|H.(j«)j - K,. , 2 2 

^ • ^(co8(2uT)+Bi.cos(wT)+B2.) +(Sin(2wT)+Bj,;^^ 



(22) 



This magnitude function ia the same for both the Butterworth and the 

Chebychev filtejE's where .i varies from 1 to N. . 

C, Chebychev Low-PaSs Fiitet 

The advantage of the Chebychev low-pass filter over the 

Butterworth low^pass filter is that the transition band of the 

respbnse at frequencies greater than w . is sharper for the Chebychev 

k c ^ 

low-pass f Uteri - This is achieved by specifying a small percentage 
of ripple in the low-pass region. The aaplitude. of the ripple is 
specified ^>y the quantity 6 (labeled RIP in the program). Figures 6, 
7 and 8 illustrate the rippling for siScond, fourth, and sixth order fil 
ters, respectively. The poles of the filter are found on an ellipse 



described by two Butten^drt^ circles of radii A and B with A<B. 
>The location of the poles on £be ellipse is a f unction of the ripple » 
6, and. ip given by the* following e<iuation: - ' ' 

VA-i((J?^-h.V/^'' t (N/e^+e'V"''^^ (23) . 

where * . \ ' - 

B Is given for the plus sign and A for the minus sign. The Chebychev 
ellipse then has major axis 3 and minor axis A. The location o£ the 
§ pl^e poles on the ellipse is given by 



I 



Real Part « A cosB , * 
Imaginary Part - B sine 



(25) 



The '6^3 are the same as given for the corresponding order Batterworth 

filter. An example of Chebychev pole locations la Illustrated In 

Figure 2. For A • 1/2 'and B - 1 in a fourth order filter, 

0 - 22.5° and 67.5"*. The, Chebychev pole legations are determined 

ftom equations (23), (24), and (25)*- 

The analog second order Chebychev low-pass filter is 

< K,a . 

H(S) - -r ^— ■ (26) 

S +KgS+K2 ' * . 



where 



1/N 



(27) 



E is calculated from equation (24) and N is the number of second 
order sections. Kg and are calculated by 



K- - 2Acps8 . ; ' (2a) 

TtKB substitution of the low-pass to low-pasa transformation for some . 

cutoff frequency u . equation (17). into equation ^(26) yields 

2 ' 

H(S) - ^, ^ ° ; . \ / (30). 

S + SKgW^ .+ 0)^ K2 
Using the extended bilinear JL transform, equation (2). and substJtutinp 
jroc for 0) we have for any section . " 

K aWDc2(Z^+2Z+l) 

«» — fc-i- ■ — — — = — . (31) ■ 

-|(z2-2Z+l) +^WDC(z2-l) + WdA2(Z>2Z+1) 

til ^ 

Qollectlng terms yields the following for the i section 

•„^,,,-.w!^ ■ (3.) 



where 



- ~ + I WDC-Kg + VfDC^K^ " . * 

aK^WDC^ \ 

hi-^ ■ ■ '''' 

2WDC^K - -f * • 



11 G, 



i varies from 1 to N. These coefficients are used tq find- |H^(j«)| 
given in (22). 

^or applications where a sharper roll-off i^ required the .Chebychev 
filters are used. The roll-off increases witV n for any fixed e. 
For fixed n, the roil-off decx;ea8es as e 'decreases. For small 
e the ripple width, 6, is small, see equation (23), but so is the 
roll-off. For larger e the roll-off improves but the ripple width 
increases. In the first case the filter will be .good at I>C and 
low frequencies, unsatisfactory at high frequencies. The converse' is 
true in the second t:ase> \ ^ 

The above observations suggest the p^rocedure to be used in 
selecting a Chebychev filter to match a set of specifications. The 
permissible ripple width sfjecifles e. With e fixed, select n 
to attain the required^oll'^f f . . v 



II. Using the- Program 

The first data card read into the program contains the number oi 
second order sections to be cascaded, N, and the type of • filter | 
desired, KN. N is equal to 1,2^ or 3, which Qorresponds to 
the 2"*^, 4^^, or 6^^ order filter respeefcively. jOI = 1 yields a 
Butterworth filter, while KN - 2^ yields a Chebychev ^lter. The 
format on the N, KN card is 212, The second data card read in * 

the sampling interval T 4n flCS format. When choosing T, 1/T 

' " * -J • • ■ 

' should be approximately equal/to ten times the cutoff frequency, o)^. 

The third data card contains the value of oj^ in FlO.4 format.. 

For the Butterworth low-pass filter, 'Is the -3db cutoff freijuency. 

" ' ' 4. 2 1/2 

For the Chebychev f ilter~ the magnitude of the response is 1/ (l+€ ) 

-1-6 at u)'» oj . U3 Is in radians. 6 is tjifi nipple factor. 

c 

If the desired filter is Chebychev, i.e., ICN -.2, the next 

data card the ripple factor (RIP) in F5.3 format. The filter 

response for all even ord^r Chebychev low-pass filters passes through 

l/n+e^)^''^ = 1-6 for u = 0 and u) . For odd order 'filters, the 

2 1/2 

magnitude Is 1 ^or oj 6 and l/(l+€ ) = 1 - 6 for w - w^. 

This program produces only even order filters. If the desired 

filter is Butterworth, i.e., KN •= ly^ tliis data card is omitted from 

* 

the data deck. 

V ■ s . 

The final data card is the starting frequency (FREQl) and the 

frequency increments (BELT) in radians. The format of the' FREQl, 

DELT card, is ZFKJ. 4. ' Determine DELT -by the f ol lowing : • 

final frequency - starting frequency 
DELT - ^Q24 

1 

This is necessary because there are 1024 frequency data points calculated 
in the program. Choose FREQl and DELT to Insure that calculated values 



will include tKa data of interest. For oaxijnuni efficiency. of the 
program, DELT should be "a awltiple pf " 2 so no deciiaal to ^i^^^V 
conversion errors are Incurred. ■ . ^ ' " 

The digital filter coefficients qre computed and printed out for 
each second order section* The full/fll^ter magnitude, response, as well 
as each section laagnltude response, is printed f or fea^h of the frequency 

' , ' . • ■■ i ' . " ■ ■ - 

increments specified. When N - 1, the section laagnitude response is 
the full filter magnitude response at^d i§ only printed opce, 

* • ' • • / ■ - ' ^ 

The program may be easily modified to ^incorporate a graphics 
display of -^the magnitude response* ^ There is ^a ii«j^ent card in the 
LPASS program indicating where the graphics subroutine c^ll card should 
be inserted. ^ . . " ' 

The ptogram Is written with Input obtained via device 4 and output 

written to device 6. These numbers should be assigned to the appro- 

*■ . * . ' . /■ 

prlate devices prior to running the program. . « 

The program was developed on the^PDP-ir/20 with a DOS/BA^CH 

operating system, frial runs" frequently uSed a TTY terminal "^s 

.well as a card reader for input (device 4) ; and a TTY terminal as 

1,-1 
' t 

well as a. line printer for output (device 6)* 

Double precision arltlmietic is employed. To decrease -required 

memory storage, only the frequency Interval values and the full 

♦ ," . 

magnitude response are saved. The section magnitude responses are 

printed; out, but are not stored. The program will produce approxi- 
mately '21 pages of output. ' 

Shown below are sample deck set-ups for the 'Chebychev and 
Butterworth low-pass filters. ^ ^ 



Data. 
Card 

1 
2 
3 
4 
5 



2 
3 
4 



Fonaat ^ 

^12 ' . 
Fib. 6 7 
F10.4,. 

F5.3 • i 
2F10.4 \^ 

212 
F10.6 
F10.4 
2F10^4 



Example * , 

0302 (3 sectioa^ Chebychev. lov-pass) 

0.001 (T - 0.001) * . • ^ ^ 

100 (u) * 100 radians) 

c . \ % 

0.10 (Ripple aiaplitude " 0.10) 

70 0.0^ (Start at u),- 70; Steps of 
Q-.06 radiai\p. Will finish just 
/past ta ■ 131 radians.) ' 

020^1 (TaectionS, 4 order, ButterwortH) 

m 

0.005 (T - 0.005) 
20 (w^ » 20 radians) 

0.0.^4 (Statt at U3 « 0, finls>' ,ni'=«f , 
past w » 40 radians in steps of 
0.04 radians.) 



-The following p^ges contain annotated examples of output data.' 

th ' 

This is an example of the output for a 4 ofder Butterworth 
low-pass filter with T- 0.0Q5 and ^w " 20 radians. The starting ' 
fr^quenpy is 0 radians and the frequency increment ^is Cr.04 radian. 



WDC = 2a. 01668 WC «= 20.00000 T 0.50000E-02 



FOR I » 1 



FOR I 



O.IOOOOOOOE+Ol 
O.IOOOOOOOE+OL. 



- 0.20000000E+01 
«0. 228697996-02 



- 0.18219614E-K)1 
2 -V^Q^- O.IOOOOQOOE+Ol 



0.83110937E+00 



B, 



O.lGOOOOOOE+01 
0.19167786E+01 



- O.2OO0OO0OE+O1 

- 0.24p59972E-02 
B2 - 0.9264O257E+00 



W 

O.ODOO 
0.0400 
^0.0800 
0.1200 
0.1600 
0. 2000 
0. 2400 



If 

O.iOGOibE+01 
•O.lOOOOE+01 
0.99999E+00 
0.10000E401 
O.lOOOOE+01 
O.lOOOOEfOl 

o.1ooooe4oi 



HI 

0.1{)000E+01 

0.10600E+01 

0.99999E+00 
d. 99997 E+OO 
0.99995E-K)0 
v0.99993E+00 
■ 0.99990E-H)0 



0. 
0. 
0. 
0. 



0/ 



H2 

10qp0E-K)l 
lOOOOE+Oi 
10000E-H)1 
l^OOOE+Ol 
OOOOE+01 

:noooE-H)i 

lOOOlE+01 



\ 



I la the stage. I varies from 1 to N. _ ^ f 

HDC is the prewarped cutoff frequency* . 

WC is the cutoff frequency. • ^ / 

T is the saffipling interval* 

X^, Aj^,.>nd A2 are the low-pass filter nuineratt)r coef f iciients. , 
B, and 9« are the low-pass filter deaomlnator coefficients. . 

is the gaiti faotor..^.. * v - - / '> -1 ^ 

W is the frequency. 

H is the overall icagmltude of the digital transfer- function. 

\, * * 

st 

HI is the magnitude of the digital transfer function (1 stage). 

nd 

H2" is the magnitude of the digital transfer functipn (2 stag§)i 

* ■■. , _ ., . ' 

H - Hl*H2. 

See Figure 4. ^ 

* ^ th ^ 

This is an example of the output for a 6 order Chebychev 

low-pass filter (three second order stages cascaded) with.T - 0.005 

' ' . ■ ■ 

and 05 - 20 radians. The starting frequency* is 0 and the frequency 

c • ■ ' . ■ ■ 

increment Is 0.04 radian.- The ripple is equal to 0.100. 

WDC - 20.01668 WC - 20,00000 T - 0.50000E-02 . ' 

A - 0.24783947 B - 1.03025433 K,, - 0.12829114 ^ - 0.99443709 
A - 0.24783947 B - 1.03025453 - 0.35049793 K - 0.5614.2438 

A - 0.24783947 » - 1.03025453 Kg - 0.47878908 K2 - 0.12841170 

FOR I - 1 Aq'« O.lOOOOOOOE+01 A^ - 0. 20000000E-K)1 ' , 

A^ - O.lOObOOOOE+pl K-^ " 0.23830688E-02 
^ -0.19774006E+O1 B^ » 0. 98727357E+00 

WDC2 - 0.40066761E+03 G(I) - 0. 16142 56^ £+06 A - 0. 96548939E+00 

FOR I - 2- Aq - OilOOOOOOOE+01 A^ - O.20d00O00E-K)l , • 

^2 " O.lObOOOOOEfOl 0. 13321469E-02 

B, - -0.1 960054 lE+01 - 0. 96557320E+00 



WDC2 -* 0.400667SiE+O3. G(I) - 0.16 30312 7 E+06 A = 0.96548939E+00 
FOR 1-3 



.*0 



WDC2 



o.oooo' 

0.040'0 
0.0800 
O.J.200 
0.160O 
9. 2000 
0.2400 



0.40066761B+O3 
H 

0.89998E+00 
0.89999E+00 
0. 90001E+00 
0. 900O9E+<m 
0v90O16E-K)O 
0.90028E-K)0 
*0.90042E+00 



O.lOOOOOOOE+01 
-O.lOOOOOOOE+01 
-0.19519^13E+01 



- o;2OO00O00E4Ol 
Rj^ 0.3O31(3788E-03 
0.95321709E+00 



G(I) • 0. 163884 96E+06 .A - 0:96548?39E+00 



HI. 

0. 96549E+00 

0.96549E+00. 

0.96530E+QO 

0.96552E+00 

0.96555E+00 • 

0.9655SE+00 ' 

0.96563E+00 



' H2 

0..96549£+0O 

0.96349E400' 

0.96551E+00 

0.96554E+00 

0.96558E+00 

0.96564E400 

0.96571E+00 



H3 

6,96547E+O0 
0-96548E+O0 
0,965476+00 
0.96550E+00 
0.96551E+go 
0.96555E+00 
X).96559E+Ou 



WDC Is the prewarped cutoff frequency. 

' *. ■ 

WG is the cutoff frequency. 

• # 

T is the sampling interval. 

B, A « Y ((^e +l+€ ) . ±(yjt +l+e ) ) 
.1 is the i stage, I varies froxp 1 to N. 



;f 



'8 



2Acos6. 



K2 « A^cos^e-fB^sin^e, /' ' 

A-., A^ , and A^ are the low-pass filter numerator coefficients, 

and B2 are the low-pass filter denominator coefficients. 
♦ 

Ik^ is the gain factor. . 
WDC2 - (WDC)^. 

G(I) = + |WDC-Kg + (WDC)^K2. 



The A following G{I) Is a 



11/2N 



Ll + e J 



W is the frequency. ^ . 

H is the overall magnitude of the digital transfer function. 



HI is tl^e nagnitude of the digital transfer function (1 stage) 

nd 

a2 Is the aiagnltude of the digital transfer function (2 stage) 

rd 

H3 Is the magnitude 6f the digital transfer func^tion (3 stage) 

• ♦ ' » 

H - H1*H2*H3. . 
See Figure 8. . ^ . ' 
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KACJSITUDE VS fIeQUENCY 
FOR f 

DIHTAL TRANSFEI^^ FUNCTION 



••• V I : 



1.0 



9^9 ' 



Oi8 



■, 3-" • • ' 
■ ■ > 

»^ 0.6 

5 '0;5 

a.2 ' 



0^707 



0' 



Jl_. — 



10 
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MA NITUDfi V8 BUEQUENCY 
buiHAL TRAKiFSi FUlJCllON 



6^" dSfc$4 BUtXERWORTH L0W-PAB8 FILTER. 




Bwt; a% U) ■ 0 fadlan 

St^pa of 0.,04 pudlaa , 
It t 3 



-4- 

20 



30 



— T— 

40 



Figure 5. 



• ^' FOR . \ 

^^t^l^AL miHSrai yUNCTION 

2HD CHSB3CCHEV LOW-PASS FILTER 



•1/(1 ♦'6^;= 1*4 



u>Q s 20 radians 
"St&rt at <d : 0 r&dian 
Eipple s 6 s 0.XO6 . 
T i 0,005 



10 



20 



(radl.ins) , 



— r- 
30 



40 



Figure 6 
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.... ■■ FOR ■ 

DIdXfAL THAN3I^ FONCTICm 



G^DfiR CHEBYCHSV LOW - PASS FXLTSR 



u»Q s 20 radians 
Sw*t at t^jt 0 radian 
Hippie 0<t 100. 

$ S D.005 

3%&pa of 0,04 radian 

-.N » ■ .... ^■ 



... " ■ :> 



—T 

30 



40 



Figure 7 
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DXSITAL TRANSFSiR FUlfCTION 



OHDSR CH£Bi6H£V LOW-PASS FILTER 




■u) 



(radiann) 



Z' 20 radians 
Start -at a s d radian 
Ripple 2 6 B 0.100 
T 0.005 

Steps of 0«o4 radian 



40 



• Figure 8 
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IHTBODPCTION . ' 

Thl8 report contains the docuaentatlon for the BPASS prpgr«ai< 
consists o£ the design propedure used, a description of the prograoi. 

and design examples using the prbgraa- 

" ■• ' ' • _ ' 

The pux^o^e of the BPASS program is" the ifeslgn of either a 

, . ■ ■ : ■ • _ : ■ -s^-. 

oaa^liaaHV flat Butterworth or a Ch^bychev filter wttlj.eflual ripple* in 
the pass band, /.^or each type of filter there is a choice of battd-pass" 
on band-stop filters . Starting with an analog filter , the bili?fcear 
Z tranBfqrm is ti^ed to desi^ an equivalent digital filter^ , Tl» user 

-ente'i-s the low-pass filter e»4er, the, type' o#>-f liter desired, then:^ .. 
sampling interval, the upper sad lower cutoff fjrequencies, the 
starting frequency and frequency increaent, and if a Chebychev filter 
Is being designed, the ripple. The low-pass filter sections are 
transformed to sifteond order band-pass or band-stop sections. Then the 
program generates the digital. filter coefflciiaits for up to six 
second order sections- in cascade or up to S 12th ordei; filter. The 

, design, is carried out in ishe frequency domain. The program calculates 
the transfer function coefficients for each second order section, the 
magnitude function for each section, and the final cascaded filter 
magnitude response over the frequency inte^al specified by the input. 

The BPASS program, written in Fortran IV is supplied as a card 
deck with this report. ' The program is in the form o^ a subroutine and 
canbe used as is by a call statement from the main program. Data 
may be input via cards with output available through a line printer. 
The Input /output devices may be altered as explained in.,this report. 
Graphic rout-ines may easily be appended to the program. 



A. PrdXlaiin^ry Ddscusslon - v 

Oae consaon b&thod of designing a digital filter is to start 
with an analog traas£d:^*£uactioa E(S) and tranafom i€ to the ^ 
digital transfer function fi(Z>. 

The transfer function, of. a second order digii^il ^Iter in the 
2 domain is giv^ti by /.^ , ^ ^ 

H(Z) » ^ o ^ ■ • (1) 

where the A*s and B^s are the coefficients of the nuo^rator and ~ 
denominator respect iVeiy . . This program will calculate ac^JL^ 
factor and the coefficients Aq, A^, Aj, B^^, and th6 

transfottaation used is the extended bilinear Z transform 

.. . • i . ' . 

wh^e _T is the sampling interval. Whetl the extended bilinear Z * 
transform ip employed, the desired frequencies must fir^t be pre^ 
warped to m^ke. them compatible with the digital filter. In the 
band'-pa^ and band-stop filterrs, the upper md lower cutoff frequencies 
/and the center frequency of the filter ate pfe interest. Calling the 
Upper and lower frequencies and oj^^ respect ivl^y, the pre* 

warped upper (WDUO , lower (WDL) , and center (WDM) frequencies and the 
bandwidth between VW and WDL are found by jfe^ 

■ ■ ... - ■■ ■ ■ I 



WB • WDU -» WDL 

■ - - » ■ • 

05 and u)v are spticif led by. thft designer and the prewairping is done 
fey the- program. 

In the design procedure for ail bandrpase W band^atpp flltete of 

^r4er n\ (n'' even), the pri^ram b^gliw by first finding ^he polea >^ 

■ ' . ^' ■ * ■ ■ , • . 

for the corresponding a* /2 order ,low-pass f liter . fSe low-paaa , 

filter Is then transformed into a band-pass or band-stop filter of 

!■ ■ ' , ■ • " .- ■ 

.... ■ 1 ,\ . 

order " *w (n' » 2n). , 

B»" But terwbrth Band-Pass Filter ' 

■ We .start with a norsaaliaed second order low-pei» ' Buttl^rWortb 

., . . ... -^^^ ^. ^ , • ■ ■ ■ ■ „ . 

filter transfer function in the S plane 

■' ■. . ' ■ ■ 

H(S) --5-^ ^ ' ■ ' . . <*> 

: , ^ s + ascpse + 1 . 

vhe re the. angle' 8 is in degrees (in the program) and isay be found 

from the Butterworth circle and the relationship 

* '- " •' , 

S - e*^^^^"' " ^^^^"^ . ' .(5) 

where n is the order of the Jow-pass flltet and "to - 1, 2, . . . , n. 
This relationahit> is determined by the followltig procedure. By 
def inltioiiV a fiUer.l?! nth order Butterworth low-pass if ^it^s 8a^n . 
chairasCteristic Is / 



1 I 2 



a 



1 + (^) 



2n 



(6) 



where a is the DC gain, is the desiixd cutoff fte^&£|E and n 

is the otder of the low-pass filter. , ^ 

."• ■ ' ' . . • ■ -. ■ " " ■ ^ . 

In the design, the poles of B(S) must be found. tl» pnJcedurfi 
is -as follows: 



|H^(jw) 



-'[H(S>H(-Sn. 



a 



2n 



2ll 



c Jw." J ' ' c 



1 + 



.2 1 



n 



0) 



.c - 



1 + 



a) 

c - 



a 



1 - 



.2 



^2 
L c J 



^' fpr .11 even 



for n odd 



(7) 



Setting the dendminators equal torero. 



(8) 



Thus, the pole locationls are the 2n roots Of ±1, depending on 
whether the loV-pasf filter order is odd or even. These roots a^e 
located on a circle with radius u)^: cefitered at the origin of the % 
plane ai^^have syoiaietry with respect to both real -and liBsginary axes. 



^0 



For n odd, a pair of roots are on the real axis and the rest are 
separated by tr/ti radians. *For n evesn, a pair of roots are located 
iT/2n fadians from the real axis and the rest Ire again separated by 
ir/n* radians. Nd roots are on the Imaginary axis, for eitt^r ev«i or 
odd n . • 



Let P]^»"«»p2n be th^Toots. Frffln the symaietry of the pole 

locations, if p3^»-«*»Pjj are the roots lying in th6 right-half plane, 

■ ■■• ' ♦ ' 

the left-half plane roots are 'f'Pn* The magnitude-squared 

functioa can then be written as . - 

2 . , vn Zn , 
a (-1) (D^ 

.H^(S)H^(-S) - (s + p^)...(S + P^)(S - Pi)...(S - p^) • 

to 1» stable, H^^S) mixat have all its poles in the left-hand plane, 
thus 

n , ' 

a<i} • 



The program is written with unity gain at DC, (u ■ 0>, therefore* 
a - 1. . 

« 

In order to locate the poles as specified above, consider the 

foJ^lowing set of equations, 

i , „e-^"^^ " , m - 1, 2,..., n; for n even 

' ' ^^^^ 

-1 - , k » 0, 1,. , n; for n odd 



Substituting equations (11) into equations (8) yields 

[L.i - .eiJ^(2tn - l)/2n , 5, . i 2..., n; for n 6ven 

- ■ - -(12) 

[S^i „ ^e*i^k/n , k - 0, l,*..-, n; for n odd 

"-u) ±k 
c 



Equations (12) will give t^e pole locations as described above » 

i 

Consider the form of equations (12) 

S - e*^^ « 0, t-cos8 ± jdlnO] . ' (13) 

FrcKQ this relationship » it can been seen that the magnitude £or each 
p^le is 0)^^ regardless of the angle, and .thus all the poles lie on a 

circle with* radius w . / ' 

c ^ 

As an example t consider a second order filter » n " 2. 
. .±i^O^ - DM . ., „ . 1, 2 

0) XID 
C 



±1 c^* 

8 - 45** . . , 

' s. ■ ^ ... 

The. relationship of these roots about the circle of radius w is 

. c ■ 

* Illustrated in Figure 1. ^ The angle 9 is always, sa^asured frrai the 

jiegative real axis. 

t^e program, only the angle (s) less than 90^ are considered 
so that the poles lie in the left-half plane because poles in the 
left-half plane are stable. Putting 8 - 45**' into equation, (4) 
yields poles at -0.707 ±J0. 707. These locations are in the left-half 
plane. From equations (12), for low-pass filter orders n ■> 1, .2,..., 
6, the values of 9 are given below.. 



Second Band-FaBS 

Loi^Pass Order \ Band-Stop 

Filter Cascaded ^ ftifi&v 

Order Angle' . Sectiona Order 

n e N n* 



1 0* 1 2 , , 

2 45* 2. '4 " 

3 60% 0* 3 6 

4 22.5*, 67*5* 4 . 8 

5 72*, 36*, 0* 5 10 

6 75*, 45*, 15* 6 12 

a is the order of the low~paas filter and is used to detenslne pole 
locations., n is also the number of second order band-pass or4band- 
stop sections which results from the transformation of the lov-paas 
filter sections mid which will be cascaded to form the band-pass or 
band-stop filters of order n*. The transformation is explained below. 
The calciU.ated angles are Incorporated in the program In the order 
given above, . - 

GlVen the nonaalleed second order l%w-paas transfer function 
equation (4)i we transform this low-*pas0 into a bmid-^asa transfer 
function for some bandwidth MB, and center frequency ^ WH by 
usin^ the transform 

SWB. 



Equation (4) then transforms to a 4th order transfer function 



H(S) 



+ 2WBCOS9 + S + WB^) + S2WB mM^ cose -1- WDM^ 

(15) 



Using the root finding subroutine '•POLRT" from the IBM Scientific 
Subroutine Package (SSP) , the roots of the denominator of equation (15) 



are found. (Note:. POLRT has been attached to BPASS as a double 
precision subroutine and is included in the card deck). The roots 
found will be complex conjugate pairs. Calling the real and imaginary 
parts of the pairs RE^^, AIM^. RE2,, AIM2 equation (15) is factored 
to yield two cascaded second order sections 

SWB SWB f^fLx 

H(S) - -5 2~ * 2 2 ' ^^^^ 

S - 2SRE^ + RE^ + AIM^ S'' - aSRE^ + BEj + AIM^ 

For each 6 of a given N, the program calculates roots for both 
sections of equation (16) and labels them the ith and the ith -f 1 
section. If N, the number of second order sections specified, is 
even, the program will calculate N pairs of RE and AIM values 
or 2N - nV roots. If N is odd, the last value of 9 is 0. 
Substituting 8-0 into equation (4) and factoring yields two 
identical first order sections, 1/(S + 1). The program will calculate 
N + 1 pairs of RE and AIM values, but because the last two pairs 
are the same dye to the identical first order sections, the last pafr 
will not be used. j 

Because both second order sections of equation (16) are of the 
same format, we will deal with only one section, the ith section 
and let 



-2RE, - D, 

(17) 

RE^ + AIM^ - 



The design of an n'th order band-pass or band-stop filter leads to 
nV2 second order sections. Substituting equations (17) into one 



section of equation (16) yields the transfer fimction f6r»the ith 
section 

h;(s) - , ^ , (18) 

S + SD^ + 

1 

The extended bilinear Z transfoiia, equation (2), is used to get to 
the digit£^ ddipain. ^ploying equation (2) on equation (18) yields 
H^(Z) for the ith second order section. 

H^(Z) - \ ^ 2D, 

Z^i^ + + C.) + Z(2C. - + C^) 

Putting the denominator of equation (19) in monic form yields the 
tranter function for the ith second order stage of the filter 

«■ 

K.,(A„Z^ + A.Z + A ) 

H,(Z)«^V5_^ ~ <20) 

^ Z + Bj^^Z + 

This equation is the same as equation (1) with the exception of the 
subsv-.-ipts. For all four filter types discussed here, the scale 
factor, Kj^, and coefficients and B2 are a function of the 

section calculated, while the coefficients A^, A^, and A2 are the 
same for all sections calculated. In going from equation (19) to 
equation (20) we have ' 



= 0 

4 2D 

T 



(21) 



li G, 



2C 



B 



8_ 

i ' ^2 



11 



B 



4 



2i 



Letting Z " e^^ « e^J'^'^ for S • joj and taking the magnitude of 



H^Cjw) we 4iave 



y^AQC08(2uT) + A^cos(a)T) + A2)^ + (AQ8ln(2o)T) + A^sln(a)T))^ 



/fco8(2ujT) + 



B, .costwT) + B2^)^ + (sin(2u>t) + B^^8in(a)T))^ 
. ^ , (22) 



The magnitude function, ^quatlon^ (22) Is ti^ same for all the filters 
discussed In this report. I 

C. Butterworth Band-Stop Filter ... 
. The design procedure la almost exactly the same as that 
of the Butterworth band-pass filter, except that the transformation 
to band-stop Is the reciprocal of equation (14), l.e. 



SWB' 



2 2 
S + WDM 



(23) 



and we find H. (S)*tobe 

After einpldying the extended bilinear Z transform, equation (2), 
we have . ' ^ 

Art - A, ^ + WDM^ 

T • (25) 

A - zmm'^ - ~ ' . . 

i r, . . •. 

and Bj^^, ^21* '^i the same functions of and as in 

equation (21). These coefficients are then used in the calculation 
of Equation (22) t6 find |H^(jw)|.^ 
° D. Chebychev Band^-Pass Filter 

^he Chebychev filter ripples with equal amplitude in the 
pass-band. The amount of ripple is specified by the quantity 6 
(labeled RIP in the program). The poles of the filter are found on 
an ellipse described by two Butterworth tircles of radii A and B 
with A < B. The location of the poles on the ellipse is a fimction of 
the rlpple'^and is given by the fQllowlng equation: 

B. A . + e-b^'" * I + -(26) 



where 



(1 - 6)^ 



'-s , ^ (27) 



and • N is numerically equal to the. order of tb0 low-pass filter 
which is transformed to yield the band-pass filter. B is given 



^ ERIC 



for the plus sign and A for the minus sign. The Chebychev ellipse ^ 
then has* major axis B anji minor axis A. -^he location of the S 
plahepoles on the ellipse is given by . . - 

Real Part « A cose 
Imsiginary Part ' 3 sine (28) 

The 8 ' s are the same as given for the corre8p<mding order 
Butterworth filter. An example of Chebychev pole locations j.9 
illustrated in Figure 2. For A *- ~ knd B - 1 in a'^rth order 
filter, e « 22.5" and 67.5°. The Chebychev pole locations a^e 
determined from equations (26), (27) and (28). v ^ 

The analog second order Chebychev low-pass filter is « 

|2/N . 

- . - * (29) 



H(S) • 



+ KgS ■ 

e is calculated from equation (27) and N is equal to the order. of 
the low-pass filter which is transformed to yield the band-pass filter. 
Ko ^nd are calculated by 

= 2Acos8 . ^30) 

= A^cos^e + B^sin^e . - 

The substitution' of the low-pass to band-pass transformation, 
equation (14) into equation (29) yields 



2 2 
S WB 



1 



2 



2/N 



•"^^^ ' + S^KgWB + ^^Ivmi?- + K2WB^) + SKgWI»f^ + WDM^ 



. (32) 



After finding' thd roots of Iquation j?32) and making the substitutions 

given by equations (17) we find the ith secW order section 

* . . ■ ■ 

swk. 



S + SD^ + 



3 



1 + 1' 



l/N 



(33) 



Applying the extended bilinear 2 transform equation (2) yields an 
equation of the form of equation (20) where - 



A a 1 



^2 - "1 



2WB S 



•(3A) 



^i ." ,t * G. 

B ' and B are the same functions of and given by 

li Zx ' 

equations (21). T^lfese coefficients are then u^ed in equation (22) to 

find iH^(ju>)|. 

n. Chebychev, Band-Stop Filter 

Given equation (29) for H(S) we apply the low-pass to 
band-stop transformation equation (23) to obtain the 4th order 
transfer function 



H, (S) 



7 2 2 

(S^ + WDM ) K, 



K + S^^WB + S^WB^ + 2K„WDmS'+ SKgWwB + K2WDM^ 

2 8 . 




(35) 



N is eiiual to the order of the low-pass filter which is transformed 



to yield the ban4-'Stop filter. 



After finding the roots of equation (35) and making the 
substitutions given by equations (17>' the it^ selcond or^er ^ectl9n is 



.H^(S) - --^ 




S + SD^ + 



1 + e 



1/N 



(36) 



Applying the extended bilinear Z transform equation (2) yields an 
^quatlon of the form of equation (20) where 



An - A„ = ^ + vim^ 
0 ^ T 



2WDM^ - ^ 



(37) 



B , and are the same functions of C and given by 

equations (21). These coefficients are then used in equation (22) 
to find |K^(ja))l. . ' 



. ^! .^a. „?i-: -r-^v 



The first data card read into the program'conta;Lns the number of 
second order sections to be cascaded^ and the type of, filter 
desired 9 KN. N Is equal to 1».2,...» or 6» which corresfK^nds to 
the order of the low-^pass filter, and hence corresponds to the 2nd, 
4th,..., oi?> 12*th order band pass or band stop filter respectively • 
KN is the type of filter desired. The values of KN specifies 
one of the fotir choices givai by 

* ^ Type y ' ' 

1 Butterworth Band-Pass 

2 Butterworth Band-Stop 

3 Ghebychev Band-Pass 

4' Chebychev Band-Stop ^ 

The foipat on the N, KN card is 212. 

The second data card -read in is the sampling Interval T in 
F10*6 fortaat. When choosing T, 1/T should be approxiiAately equal 
to ten times the center frequency (WDM) • 

The third data card read in contains the values of the upper 

an,d lower cutoff frequencies, u) and u)- , In 2F10.4 format. For 

• u 1 . - 

the Butterworth filters, the cutoff frequencies are the -Sdb cutoff 

frequencies. For the Chebychev filters, the magnlttide of the 

2 ^ 

response is 1/(1 + e ) » 1 - 5 at the cutoff frequencies, u) is in 
radians, is the ripple factor. 

If the desired filter Is Chebychev, I.e., KN » 3 or 4, the 
next data card coi>talns the ripple (RIP) factor in F5-3 format- 
If the desired filter Is Butterworth, i.e., KN = 1 or 2, this 

« 

card is omitted from the data deck. 



The final d^ta card is the starting frequency (FREQl) and the 
frequency increments (DELT) in radians. The format of the FREQl, 
DELT card is 2F10.A. Determine DELT by the following: 



PELT - final firequency|^-' starting frequency 

This is necessary because there aVe 1024 frequency data points 
calculated in the ptogram* Choose^ FREQl and DEI»T to Insure that 
calculated values will include the data of 'interest. For maximum 

— K 

efficiency of the prpgram, DEl-T should be a multiple of 2 so no , 
decimal to binary conversion errors are -incurred. 

^ The digital filter coefficients are computed and printed out for 
each second order section. The full filter magnitude response, as we3 
as each section magnitude response, is printed for each of the 
specified frequency increiMnts. When there is only one second order 
section, the se(^^)n magnitude response is- the full filter magnitude 
response and is only printed once. « 

The program may be easily modified to incorporate a graphics 
display of the magnitude response. There is a comment card in the 
BPASS program indicating where' the graphics 8ii}>routine call card 
should be inserted. 

The program is written with input obtained via device 4 and 
output written to device 6. these numbers should be assigned to the 
appropriate devices, prior to running the. program. 

The program was developed on a PDP-11/20 with a DOS/BATCH 
operating system. Trial runs frequently used a TTY terminal* as well 
as a card reader for input (device A); and a TTY terminal as well 
as a line printer for output (device 6). Double precision aritrhmetlc 



la esployed. To decrease required memory storage, only the frequency 
LntJrvnl vn^ues and the full magnlttide response are saved. The 
section fflagnittide responses are printed out, but are not stored. 
The program will prbduce approximately ?1 pages of output. 
Shown below are gasiple deck set-ups. 



Data 








Card 


Format 






1 


212 




0504 {5 section . Chebychev ban -h 




F10.6. 




0.002 (T - 0.002) 




2F10.A 




60 40 (u - 60, 0) • 40 radians) 

U 1 ' 


4 


F5.3 




0.10 (Ripple amplitude - 0.10) 


.5 


2F10.'4_ 




0 0.1 (Start, at w » 0. Steps of 








o!l radian. Will linlsh 








Just past to - 102 radiar- . 


1 


212 




0401 (4 sections Butterworth hai.f1 




Fm.6 




0.002 (T - 0.002) 


3 


2F10.4 




60 40 •= 60, a)^ - 40 radians) 


4 


2F10.4 




0 0.1 (Start at w - 0. Steps of 



0.1 radian. Will finish 
Just past (D - 102 radians). 

The following pages contain annotated examples of output data. 



This is an exain|^e of the output for an 8th order Buttetworth 
band-stop filter (-0402) with T - 0.002, w - 60 radians, aad 

• ■ ' ■ ■ 

(0 » 40 radians. tJie starting frequency is 0 radian and the frequency 
increment is 0.1 Vadian. 



WKJ « 60.07210 HDL = 
T - 0.20OO0E-O2 



40.02135 mi 



49.02902 WB - 20.05076 



THE ROOTS OF THE FILTER ARE GIVEN BELOW 
REAI 1) « -8.52659476 IMAGINARY(l) « -44.46786740 
RL.L I) - -9.9^/88916 mAGINARy(2) " -52.140959^; 
REAL(3) « -4.S3076557 £MAGINARY<3) - -59.01588870 
REAL 4) - -3.12232691 IMAGINARY(4) - -40.49140500 



FOP 



lEFFICTENTS OF EACH DIGITAL FILTER SECOND ORDER SECTION ARE 
GIVEN BELOW . ,/ 



- 1 A-^ = 0.10024038E+O7 A. 
A^ = 0. 100240 38E+07 
- b: =f -0.19584863E+01 B, 



-0.19951923^+07 
0.9812548iE-06 
O.96653295E+00 



FOR 



2 A 



0 



B, 



0.10024038E+07 
0.1O024038E+O7 
-6.19498774E+01 



b; 



- -0.19951923E+07 
« 0.9t769448E-06 
« O.96090O48E+OO 



FOR 1 



FOR . 



W 

O.OC )0- 

0.1000 

0.2000 

0.3000 

0.4000 

0.5000 



3 A„ = 



0 



B 



1 



H 



0.1OO24O38EH-O7 A 

0.10024038E+07 K: 

-0.19681836E+01 B2 

0.10024038E+07 A 
0.10024038E+07 

0.1',\ Z". 01 B„ 



O.IOOOOE+Ol 
O.lOOOOE+01 
O.lOOOOEH-01 
0.99999E+00 
0.10000^+01 
O.lOOOOE+01 



=. -0.1995 1923Ef07 

- 0.98755180E-06 

- 0.98202 35 3E+00 

*'-0.19951923E+07 

" 0.99216788e4)6 ^ 



Hi 

0.11726E+01 
0.11726E+0r 
0.11726EH«1 
0.11726E+01 
0.11726E+01 
0.11726E+01 



"f 



\ H2 
0.a5284E+(K) 
0. 85284 E+00 
0.852 84 E+00 
0.8S28^E+00 
0.85283E+00 
0.852 82 E+00 



H3 

0.68611E+00 
0.68611E+0Q 
0.686UE+00 
0.68610E+00 
0.68609E+00 
0.68609E+00 



H4 

0.14575E+01 
0;14575E+01 
0.14575E+O1 
0.14575E+01 
0.14576E+01 
0.14576E+01 



WDU is the prewar ped upper frequency* 
WDL \b' the prewarped loweV frequency. 
yM! is the prawarpe'd center frequency. 
WB is. tli^ bandwidth, WDU - WDL. 

T' is the sampling interval. . . - 

The Real- and Imaginary part of the roots of tlie filter are ^iven next. 
I is the ith stage. I varies from 1 to N. , ^ 

A^, AJ . are the Butterworth band--stop filter numerator coef f icientSw 

012 ^ - * 

is the galn^ factor . " ^ 

B , and are the Butterworth ban4-8^op filter denominator coefficients 

1 o 2 • ' 

W is thd^ frequency - " • . 

H is the overall magnitude of the digital transfer function 

HI. is ^he magnitude of the digital transfer function (1st stage). 

H2 is the magnitude of the digital transfer function ^(2tlff^a^ge). 

H3 is the magnitude of 'the digital transfer function (^rd sta^e). 

H4 is tfe magnitude^ of the digital transfer function (-^h stagfe). 

• * ^ ' ^ ^ ' ' 

■ i. 

See Figure' 4. 



Th^s is an example of the output for an 8th order Ghebychev 
band-pass filter (0403) with T = 0.002, w = 60 radiana, and 

• * * 

0),'* 40 radians. The starting frequency is 0 radian, the frequency 

1 . . • , . ■' • - 

»■..■.-. . 

increment is 0.1 radian, and . the ripple is 0^,1. 



WDU i= 60/07210 WDL » 40.02135 WDM « 49.02902 WB « 20.05076 
= 0.20000E-02 , 



'a = 0.37642105 B =%1. 06850027 K.^ = 0.69553541 IC » 0.28813942 
A = 0.37642105 B - 1,06850027 Kg = 0.28810020 - 0.99524620 



THP ROOTS OF THE FILTER ARE GIVEN BELOW 



REALCD =« -3.19528085 
REAL(2) = -3.77772492 
REAL(3)= -1,15829660 
REAL (4)" = -1,73001696 



IMAGINARY (1) 
IMAGINARY (2) 
IMAGINARY (3) 
IMAGINARY (4) 



-44,97792290 
-53.17662810 
-40.10115290 
-59,89456990 



THE COEFFICIENTS OF EACH DIGITAL FILTER SECOND ORDER SECTION ARE 



FOR 



FOR I 



FQR I = 



FOR 1 



GIVEN 


BELOW 


1 = 1 


^2 = 
h = 


I - ? 




1=3 


h ^ 


1=4 


A = 

? = 

2 

^1 = 



0,10000000E+01 
>-0. lOOOOOOOE+01 
-0,19737935E+01" 

O.IOOOOOOOE+Gl 
■0,lOOOObOOE+01 K - 
^ -0.19889723EfOi Bj » 

o.iooooodOe+oi a, = 



A * 0,0QO00O00E+(X) 
= 0,10395603E-01 
B2 = 0.99732565E+00 

A. - O.OOOOOOOOEH-00 
« O,10375296E-01 
hi = 0.98504460E+(K) 

A = O.OOOOOOOOEf(» 
K, - 0,19406846E-01 
0»9953'84,93B+-00 

O.OOOOOOOOEfOO 
= 0.19346636E-01 
b: = 0.99312838E+00 



• w 
6.000 
0.1000 
0.2000 

6.3000 
0.4000 
0.5000 



0 . OOOOOE+00 
0.12493E-12 
O.19990E-11 
0.10121E-10 
0.31991E-10 
O.78116E-10 



HI . 
O.OOOOOE+00 
0.51560E-03 
0,10312E-02 
0.15468E-02 
O.20625E-O2 
0.25783E-02 



H2 

O.OOOOOE+QO 
0. 36886 E-03 
0:73773E-03 
0.1 1066 E-02 
0.14755E-02 
0.18445E-02 



H3^ 
O.O.OOOOE+00 
0.12106E-02 
0.24211E-02 
0<36318E-02 
0.48426E-02 
b.60536E-02 



H4 

O.OOOOOE+00 
0.5426.5E-03 
0.1O853E-02 
0.16280E-02 

0.2l707E-^02 
0.27134E^02 



WDU is the prewarped upper frequency. 
WDL Isi the prewarped lower freijuency. 
ytlM is the prewarped center frequ^cy. 
WB Is the bandwidth, WDu!- WDL. . 
T is the sampling interval. 



n A 1//. / -2 . , , -1.1/N ^ , ^ ~1.-1/N. 

B, A • 2'^(/e + 1 + e ) ' ± (/£ + 1 + e ) ) 

Kg - 2Acos(e) • .. .. ' 

^2 - A^cos^O) + B^sin^(e) 

* , 

The Real and feiaginary part of thie roots of the filter are given next. 
I is the ith stage* I varies from 1 to 

Aq, A^^ A^^s^ the Chebychev band-pass filter numerator coefficients. 



e cairN 



is the gainAfactor. 
B^, and B* are the Chebychev, band-^pass filter denc»alpatpr coefficients. 
W is the frequency, 

H is the overall magnitude of the digital transfer function. 

HI is the magnitude of the^digltal t^ransfer f uric t ion (1st stage). 

H2 is the Ti^gnitude of the /i:^ital transfer function , (2nd stage). 

H3 Is the magnitude of the digital transfer function (3rd stage). 

H4 is the magnitude ^f the digital transfe/r function (4th stage) • 

i ^ / ' . ' ' 

See Figure 5 . * 




M^iGNITUDE VS FKEQUENCY 
FOR 

DIGITAL TRANSFER FUNCTION 



8^^ ORDER ^UTTERWORTH BAND-PASS FILTER 



0 radian 



N - 4 

Start at u) 
T « 0.002 

Steps of 0.1 radian 
u) 60 radians 



u 



0)^ 



« 40 radians 



1.0 



ad 



-a 

3 



00 

m 
?1 



.0.8 
0. 7 
0.6 
0.5 
0.4 
O./i 
0.2 

0. 1 

(/To 



0 




a) 



(radians) 



■J 9 



100 



Figure 3 



5 MAGNITUDE. VS FREQUENGY 

FOR 

DIGITAL TRANSFER FUNCTION 
8^^ ORDER BUTTERWORTH BAND-STOP FILTER* 



N « 4 . 

Start at u> « 0 radian 

T « 0.002 . 

Steps of O.I radian 

u) = 60 radians 
u 

03- « 40 radians 
X - ' 



1.0 



o 

01 
TS 
3 

•H 
C 
DC 



0.9 - 



0.8 



0. 7 



0.6 



0.5 



0.4 



a. 3 



0.2 



0. 1 



0.0 




0 



40 



60 



80 



lOO 



CO 



0) 



1 u 
a) (radians) 



Figure 4 



MAGNITUDE FS- FREQUENCY 
FOR 

DIGITAL TRANSFER FUNCTION 
8^^ ORDER CHEBYCHEV BAND-PASS FILTER 



N • 4 

Start at ui r 0 radian 
Ripple « o 0.100 
T - 0.002 

Steps of 0.1 radian 
u) » 60 radians 
uj" = 40 radians 



•ft 

5=: 



O.i- 



0.7- 



0.6- 



0. 5- 



0.2r- 



0.1- 



o.cL 




\ 



100 



fi) (radians) 



ERIC 



.6 



Flg&re ■> 



MAGNITUDE VS FREQUENCY 
FOR 

DIGITAL TRANSFER FUNCTION *' 



10 ORDER CHEBYCHEV BAND-STOP FILTER 



N - 5 

Start at 0) « 0 radian 
Ripple > a = 0.100 
T 0.0G2 

Step? of 0.1 radian 



0} 

I 

05, 



U 



60 radians 
40 radians 



o 

;;3 



to 



r 0.8^ 




(radians) 
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INTRPDUCTION 

This report describes hnw to use two programs for the weighted 
least squares desi^ga of ^aohrecurslve^.and recursive digital filters. 
First the theoretical aspects are considered and the design equations 
arl developed* The signal imidel in this work Is assim^d, to be a 

polynomial because the state model Is sltaple. Hdwever, the ^Keory Is 

■ . . ' ; , . 

easily extended to Include any signal model represented by a linear 
differential equation. - , 

Then the operation of the two prograas^ 1^ described along with 
examples to ll!|.u8trate their operation. 




LEAST SQUAKES DESIOJ OF NONRECURStVi: DIGITAL FH/TERS 

The design of nonrecurs^tve and recursive digital filters using 
weighted ^^east squares is based on a model for the Input signal. 
The most coranon models that are currently in use are^dlf ferent lal 
equations « The subsequent representation of the diff erentitrl equa- 
tion by a flrst-^rder vector differential equation leads to the 
concept of a state variable and the state space representation for 
a system. By use of the proper fort&ulatlon» a continuous signal 
represented by a differential equation can be described by a first- 
order vector difference equation or a discrete tin^ state space 
model* References that describe the essential' aspects of state 
variables are by De Russo^ Roy and Close [1] and Chen [*2]. 

Since signals from laboratory Instruments are not usually des^ 
cribed by a differential equation, approximatlonH often utilize a 
polynomial* For this reason, the vector form of a polynomial approx 
Imatlon will be used in this report. The reader should he aware 
that this can be generalized to include any signal irodel that can be 
represented as a linear time-varying differential eqitotlon* Later 
in the paper a scalar irodel representing a Gaussian time signal 
will also be utilized to develop a time-varying filter that can be 
used for reducing the base-linfe error and for the initial separation 
of signal components. * 

To develop the polynomial itodel let the signal z(t) be repre-^ 
sented by a polynomial of order mT^ At the time t =* nT the state 
vector for the signal is glyeQ by 



(1) 



Redefining the state vector as 



t"nT 



x(nT) ^ 




(2) 



t-nl 



the use of a Taylor series representation for each elenient of x(nT) 
now permits the state of system at t-(n+h)T to be described in terms 
of the 6tate at t««nT by the relationship . " 

x{(n + h)Tl - *fh"lx[nTl (3) 

s 

where hi is the mxm state transition matrix with elements 



$ fhl « iL_ 



1-1 
h-^ 0<l<m 



(4) 



The state transition matrix i\h] satisfies all of the rel^t lon«hI pn 
for general state transition matrices with the important t^belng foi 
this work V 



*[-h] - [*(h)J 



-1 



(5) 



-2- 



r o 

ERIC 



At this point, th&se laadels can be utilized in the design process* 

Design of Noprecursive Filte rs 

Let the input signal start at time t«0 and assume that the signal 

over a finite data window is approximated- by a polynomial 2(t), 

■ o 

Defining the state of the signal by (2), then the state of the signal 
at time t=»(n+h)T is given in terms of the signal at time t»nT by (3). 
Given that the signal starts at t-0, the first I observations 
are each defined as ^ 

■ I 

lUT] - Mx[iT] + vljT] j«0,i,2,...£.-l (7) 

where M is a row oatrix that relates the measurable state variables- 
to the actual measurements. The elements of each noise vector yijT]' 
are the measurement noises. In general these arfe taken as random 
variables with zero mean and the time dependent ^tbcovarlance maltrlx 

<^ RljT] - E[v(jT)v^jT)-] (8) 



however, in most laboratory 'systems "only the data ife-available and 
the deysivatives are not measurable. Furthermore* the noise covariance 
matrix is usually not known and jfor scalar measurements the noise 
variance is assumed constant and ^he noise sampl^ iin^rrelated. For 
the remainder of this paper, only scalar measurements and uncorrelated 
measurement noise with time-invariant statistics will be assumed. The 
mean value of the noise will be taken as zero and tlie variance as 0^. ^ 
The results are easily extended to vector measurements and to measure- 
menl? noise with time varying statistics. 



^'8 



For £ observations, the tota,l observation vector at t«^Tvi8 

' defined as ^ * 

% 

* -mm 

ia[nTl 
JBl(n-l)T} 



(9) 



M[(n-£+l)T|j 

This vector n<«f forms a data window of 8. data points. For an estimate 
of the data at t-nT the' use of the expression 




* y x[(n-j)T] - <&(-J) x[nTl, 

is combfned with (9) to yield 



M 

M*(-l) 



xCnTl + v^lnT] 



\ 



The matrix of constants H[nT] is now defined as 



• H[nT] 



•M_ 



s^that (11) can be' written as 



(10) 



(11) 



(12) 



m^lnt] - HfnT] x[nT]^ v^[nT] 



(13) 



Tlie elements iiji HinT] are constanta given by 



9 o 

ERIC 



H 



(-1)^ 0<i<l 
-4- 



(14) 



where Z is the size of the data window and is the order of the 
polynomial 'plus one. When \-j-0 the value* of (14) is one. Note that ' 
slxxce (12) is a laatrlx of .-constant's there Is no ii€?ed to mallei the 
matrix a function of time. However, jja the derivation foe. the recur- 
sive filters this provides a method for) separating different H matrices. 
The optimal estimate of the data/at t-nT is now given In terms 

of weighted least squares or minimum variance because this form is 

* '. ■ 

Utilized in the derivation of the recursive filters, tf the covari-^ 

I ^ ') 

ance matrix for the total observation vector is 

R^(nT) « Elv^ (nT) v^(nT)j (1$) 
the optiiaal minimum variance estimate is - 

xInT]'- W[nT] m .[nT] ' • ' (16) 

where W[nT] is a series of constant weights giv^ by 



W{nT) 



H^(nT)[R^(nT)]"4(nT) 



'V(nT)lR (nT)l~^ » (17) 



2 

For uncgrrelated noife with constant^^lance o (15) is a diagonal 
mj^trix with elements and (17) redSSes to 

^ W[nTl - [H^(nt) H(nt)]"^ H^(nT) • " ' (18) 

This is the same result obt4in^d using conventional least squares when 



the noise covarlance Is a^^^S^gonal matrix of equal constants. The 
reader should be awsiN^-that the estimate vector In an estimate of tht^ 
^^4ata and air of the derlvdtlvea In the model. However, all of the 
\y\ weights in the derivative terms must be scaled by the scale factors in 
the state vecrtor define^ by (2). 



The covariance matrix for the estimate given by (^.6) is given by 



S(nT) - W(nT) RJnT) Wt(nT) 
Stibatftutinjg (17) lotb (19) nov yields 



(19) 



S(nT) - iH^inT) {R^(nT)]'^ H(aT)]~^ 



idilch sliQ>ll£les further to 



(20) 



S(nT) - [Ht(nT). H(nT)]~^ (21) 
for the uncorrelated noise. 

By defining a delay or prediction factor a, estimates of the dAta 
aT units behind or ahead of the point t«nT can be done. To do this 
the^tal observation' vector is written as 



m^(nT) • H^(nT) x[(n-a)T] + v^(nT) 



where H^(nT) now takes the form ' 



H^(nT) 



M*(-fx) 
M$(-a-l) 



M<^(-l-a+l) 



■ L 

The individual elements of the matrix now become 



(22) 



(23) 



H^(nT)^j - (a-i)J 0<i<ll 

0<J<q 



The form for the optimal estimate nov utilizes H (nT) in (18). The 



(24) 



n 



covariance of the estimate uses H (nT) In (21). 

a 
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•For a f iv€--po:|jnt window with a«0, the optimal weight Mtrlx and 
the cqVarlance matrix for the estimate are shown in Table 1 for. a 
third-border polynomial fit. For the estimate of the data in 

the middle of the window is obtained and weight matrix and covariance 
matrices are shown in Table 2. This case corresponds to weights 
given in reference [3l« 

In practice, the design of npnrecurslve filters is initiated by 
specifying the si^e of the window which is the number of rows In the 
H matrixt the order of the polynomial approximation which Is one less 
than the number of columns in the H matrix and a which determines 
the coefficient values in the H matrix. ^Thls makes the design suit- 
able for use with interactive graphics since only three parameters 
need be specified to generate the weight and covariance matrices. 

If the model of the signal process is modified by additive 
uncorrelated driving noise ♦ the variance terms of th^ driving noise 
fade the lieraory so that past data has less ^ffect on the estimate. 
This can be thought of as uncertainty in the signal model. The con- 
cept is particularly important In recursive filter design^ For non- 
recursive filters, it can also be utilized and can be useful when 
using a non-recursive filter to Initialize a recursive filter. 

If the model includes driving noise, the state at time t<mT la 
given by 

xlnT] - 1>(-11 xI(n-l)Tf + wI(n-l)T] (2 

where W[(n-*1)T] is a sample from a noise process with mean zerc^nd. 
covariance matr^ Q. This noise process is assumed to be white. In 

-7- IHL 



terras of the to&^l measurement vector, a^JnTj now becoiQes 

^[i^T] - H[nTl x[nT] + 2^(nT) + v^[nT] 
^ere £j.InT] Is the total noise vector due to the driving noise 
scalar measurements this Is given by 



(26) 
For 



£^[nT] - 



0 

-M*[-l] wI(n-l)Tl 
-M<&(-21 wr(n-l)T] 



- m[-l] w[(n-2)Tj 



(27) 



The total noise vector that corrupts the total measurement vector is 



now defined as 



(28) 



For scalar measurements the covariance matrix for r^[nTl has diagonal 
elements whose value increases down the diagonal. Since the diagonal 
elements ar^ not the same, the minimum variance expressions of (16) and 
(17) must be employed to find optimal linear estimates. 

Example . • 

to illustrate how the driving noise Effects the filter weights, 
consider a zero order process which is equivalent to estimating a signal 
of constant value. For a three-point filter with scalar measurements, 
the total noise/vectoi- is 



f" 



r^tnTj 



V ( nT } 

v[(n-m| 

vI(n~2)T] 



wl(n-l)Tj 

w[(n-2)T] - »f (n-l)Tl 



(29) 



1 



If the variance of the s^asurement noise 1^ and of the driving 
2 

noise Oj^ the covarlance matrix of rj.InT] is •» 



R^EnTl 



0 



0 cj2 + o 2 0 



0 



2 2 
o + 20^ 



TKe resulting weight ©atrlx obtained using mlnlraum variance Is 



(30) 



where 



(o^^ + o^) (20^^ + ah 



o (2o^ -f- a^) 



» 



(31) 

(32) 
(33) 



and 



o^ff 2 4- a^) 



•2 4^ 



(34) 



^ where D la given byf 



D - 



(20j^^ + o^> (o^^ + o'^) + 0^ (20^^ + CT^) 
+ a2(a,^ + 0^) *T 



(35) 



The terms In the weight matrix satisfy the following Inequality 



^0 ^ S - ^2 



1/: 



(36) 



For 0^ «'0, the equalities hold atid * 82 1/3 which are 

2 

•the well known weights to estimate a raean. If 0^^ Is not zero, the 

2 2 
inequalities hold and as becomes large with respect to a , 

* 1* 
approaches one and a^ and a^ become smaller so that fading Is 
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\ 
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introduced. The covarldnce of the estimate is a scalar. given by 

? (o 2 + a^) (2(} 2 

S(nT) -- Co2 (37) 



where C ia a dia»asionXess constant whose value approaches one as a\y^ 
becoToea larger than a , .Thus for » a only the current observa- 
tion ia effectively used and the variance of the eaticdte is approximately 

2 : ■ ' \ ■ ■ . 

The fading obviously reduces the signal-to-noise enhanceo^nt of a 
nonrecursive digital filter. On the other hand, fading can be used to 
reduce ,t he deterministic error due to a. nonexact model. In practice, 
fading in nonrecursive filters^is not often utill^zed. Instead, the 
window length is more typically used as a desigt^ parameter. In future 
work, however, it may be desirable to further explore the relationahlp 
between fading and the aize of the data window to achieve Improved designs 
when nonexact signal nxjd^ls are employed. 

LEAST SQUARES DESIGN OF RECURSIVE DIGITAL FILTERS 
The fixed memory or nonrecursive' filter design is now extended to 
recursive digital filters that utilize all of the data. The resuU Is 
a recursive form that is usually called the Kalman filter. While there 
are a variety of derivations for the Kalman f 11 t«r, stairtlng with a fixed 
memory filter using a polynomial model with driving noise gives the result 
in such a way to give the reader greater intuition about how the filter 
works. 

i 

The derivation of the recursive filter Is started by using a signal 
model • ^ 

/ 

x[nTj - <{•[!] x{{n-l)Tj + w{(n-l)T] ^ (38) 



and the scalar observation or aieasurement i^odel 

yCnT] - MxInT] + v{nTj. • (39) 

. / 

In (38) and (39) the terms wrXn-DT) and vln^r| art- the dMvrng iioIh.. ' 
and the measurement noise. These noise terJa have zero mean and are 
uncorrelated with themselves and each other. Given the n+l^Bjeasureiaents 
starting at t-0, the total observation vector at time t»nT Is given In 
terms of x{nT] as 

XfnTi » H[nT] ^ £j.InTj + y^[nTj - (40) 

In (AO) the matrix HlnT] is 



) 



H[nTl - 



the vector j^fnT] is 



M 



1 



MH-n] 



(41) 



1 



0 . - 



-M*f-1} w[(n-l)Tj 



1 



-M^Sj *[-(3-J)) w[(n-J)Tj 



-M^?i *[-(n+l-j)] w[(n-j)Tj 



(42) 



* and v^{nTJ is the total measurement noiae vector. Defining the , sum ^ 



r^fnT] - p^[nT.l +.V^InTl 



(43) 
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in 



the, total noise vector vith autocovar lance 



\[nT] - E{r^[nTJ r^InT]] 
the optimal estimate using linear jninlmuia variance is given by 



(44) 



x[nT] - 


i 

W[nT] \ 










where WfnTj is the weight matrix 










W[nTj - 


H^[nT3 











(45) 



The covariance of the estimate is 

S[nT] - W[nT] R^[nTj W*^[nT] 
Substituti<$n of (46) and (47) gives an alternate form 



1 

(46) 
(47) 



S[nT3 



H'^lnTj 



-1 



HinT] 



which allows the alternate form for (46) to be writ 



-1 



ten as 



(48) 



W[nT3 - S[nT] H^(nT] 



[nT] 



-1 



'(49) 



the forms givA by <46) , (47). (48) and (497 apply to the remaining 
estimates used in the derivation and &jie^ reader should remember them. 

Next the prediction or forecast of the signal state X[(n+1)T] is 
found by first writing the total observation vector X^.^""^! a® 



1 



Y^lnT] - H^[nT3 x[(n+l)T] + p^^Jnl] + v^[nT] <50) 



where H- [pt] is given as 

-L 



H^(n^] - H[nTJ <^[~1J 



(51) 



and £ 



.^Jnl/j is 



R+1 X n+l matrix with elei^nts a , Taking the covariance of p7.[nTl 



-t 



and performing aoma algebraic manipulation gives / 

/ 

E{£j^^.[nTl iJjnTl) - E^InT] £j[nT]} + H^[nTl QH^tpt] (57) 

, - ' _ ; / ' 

where Q ia the diagonal covariance matrix of the driving noise vector 
v[nT]. This matrix is usually asaus&ed to be time ^^nvariant. Thus^ from 
(57) - . \ 

R^^[nT] - R^[nTl + Hj^[nT] QHj[nT] (58) 
Subst^-tutlng (58) Into (36) now gives 

Sj^t(tt+1)T] - Wj^[nTl R^InTl wJ[nT] — (59) 

+ W^lnTj H^[nTl QH[[nT] wJ[nTl 



ibiasedVstj 



If X-{(n+l)T] is an imbiased^estizoate then the Constraint relationship 

W^[nT] HjlnT] » I (60) 
must be satisfied [4] so that (59) can be simplified to 

' . S^[(ii+1)T] « Wt^tnT] Rj.InT] wJ[nTJ + Q (61) 

^Also Recognizing that x [(n+l)T] can be written as . 

X [(n+l)T] = #[1] xInT] \ (62) 

the weight matrix W^^InT] is ^ 

W^tnT] * *[1] W(nT). (^3) 

Substitution at (58) intt> (61) and app^j^lng (47) now gives the 'final 
form for S [(n^^lXT] as 

S ((n+l)Tl - «l>ll] S(nT] *^{nT] +Q (6A) 

The recursion Ih now formulated. When the obsfr%?atlon at t * (n+l)T 
arrives the new total observation vector in terms of the signal* state 
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jti (iH-l)T] is N 

Ztl("+1^T] - Hl(n+1)TJ + £^f(n+l)1') + v^{(n+l)T). 
The estimate la given by • 



x((n+l)T] - Sl(n+1)T] H^[(n+1)TJ R^Kn+DTj 



-1 



with the CO variance 

S[(n+1)T] 



H^'LCn^-DT] [R^Kri+DT}]""^ H[(n+1)T] 



-1 



(65) 



(66) 



(67) 



where R [(r4-1)T] is the covariance matrix of the total measurement noise 



vector 



r^[(n+l)T] - £^[ (n+l)T3 + v^[(n+l)lj] 

First the recursion for the covariance taatrlx 

-M«[-13 wInT] 



(68) 



is found. Given that 



£^.[(n+l)T] 



-M I *[-(3-.1)] w[(n+l-l})T] 
j-1 



n+1 

-M Z *(-(n+l-j)] w[(n+l-J)T] 
j-1 



substitution of (^2) into (69), gi 



ves 



£^[(n+l)T] 



Next HUn+l)T] is given as 



0 



(69) 



(70) 
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H[(n+1)T] - 



M 

M*[-l] 



|^M(l>I-(n+l) 1 



^'9 



(71) 



^ich can be rewritten with the aid of (41) as 



7 



Hint; 



M«(ii 

H[nT] 



(72) 



The covarlaiice matrix for v^[{n+l)Tl la now an n+2 x n+2 diagonal 

2 ■ . . ■ ^ 

matrix with elements 0 so that Rj.[(n+1)T1 Is seen to be 



R^[(n+1)T] 



a^l 0 



0 I R^^[nT] 



(73) 



Since this Is a diagonal matrix Its inverse is 

T 



R^[(n+l)T] 



2 I (T 



■ r r 



0 



-1 



(74) 



Substitution of (74) and (72) int0 (67) and performing the matrix 
niultlplicatlpn now yields 



S{(n+1)TJ 



^~ + l*^!-!] H^[nT] 

o 



•X 



H[nT] 



-1 



(75) 



Substitution of (51) ipto (75) and the use of the inverse of (56) gives 



S[(n+1)T] 



-1 



-1 



(76) 



Equation (76) along with (64) now forms a recursion for the covariance 
of the estimate. 
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SFor the quantities In (76) the application 6f the matrix Inversion 

lemma [5J, gives * 

-1 



(77) 



This form will be used to generate an exptssslpn for the Kalman gain. Post 
multiplying both sides of (75) by M^/a^ and using (75) the Kalman gain is 
defined as " ^ 



K[(n+1)T] - S[(iH-l)T] M^/o^ ' 

Xr 2 



(78) 



S^[(n+l)T]-M'-[0^ + M Sj[(n+1)TJ M*']"^ 
Thus the covariance loatrlx given by (76) becomes 



S[<n+1)T] - [I - Kt(n+I)T] Mj S^Kn+DT] (79)» 

The recursion for the optimal estimate is now formed that uses the Kalman 
gain given by (78). Using the form for the optimal estimate given by 
(47) the estimate xl(n+l)Tj Is 



x[(n+l)T] - S[(jnfl)T} H!^[(n+1)T] 



R^j.[(n+1)T] 



-1 



y^[(n+i)T]> 



(80) 



where y^[(n+l)T] is the new total observation vector given by 



y[(n+l)TT 



X^[(n+1)TJ - 
I 



Substitution of (74). (72) and (81) Into (80). and carrying out the matrl; 



(8i) 



multiplication yields 



I 



x[(n+l)T] • Sl(n+1)T] 



M_y[(;n+ijri 

2 

a 



(82) 



+ *C-1] H^In^) 



R,,rnT] 



-1 
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Using the estimate of the forecast x. [(n-fl)T] (82) Is simplified to 



M 



x[(iH-l)Tl, - SKn+l)Tl iij yl(n+l)Tl + |SjI(n+l)T} 



[J 



'I 



x^.[(n+lXTl 



(83) 



Adding and subtracting (M M/a ) Xj^l (n+1)T], performing some algebraic 

< 

manipulation, and applying (76) now yields 



x[(n+l)T] - Xjl(n+1)T1 + K((n+1)T] 



y[(n+l)T] - M x,[(n+l)Tl 



(84) 



where K[(n+1)T] is given by (78). Equations (84), (79), (78), (64) and 
(62) now forfc the recursive formulation called 'the Kalman filter. These 



equations are now summarized as 



and 



Xj^(nT) - $(1) x[(n-l)T] . 



S (nT) - *(1) S[(n-l)tl <I>^(1) ,+ Q , 

K(nT) - .Sj,(riT) H^[a^ + M S^(nT) m'^)"^ , 

^ ■ . 

S(nT) - [L-'^lInT) M] S, (nT) 



/ 



(85) 
(86)' 
(87) 
(88) 



/ x(nT) m Xj^(nT) + K(nT) [y(nT) - M Xj^(nT)] 

^ fl - K(nT)M] x[(n-l)Tl + {C(nT) y(nT) 

In this set of equations^ x, (nT) is the forecast or prediction of the* ' 

estimate at t«nT using the prevlou«ly generated esjtlmate at t«((n-l)Tj. 

The covAriarice of the forecast Is'^^CnT). The tepn K(nT) is ^he time 
*^ ' ■ ^ - 

varying Kalman gain matrix and y(nT) is the observation at t*nT^. The ' 

terms x(jj ^ and S(nT) are the estimate* at t«nT and *lts covarlance^ 



(89) 



The term a is the varianc^ of the measurement noist? and the term Q lf4 
the c«^riance of the driving noise. It is the term Q that serves 



a key design parameter. 



If there is no driving noise so that the diagonal eleroentB \)f Q 
are all zero^ the fljlter Is simply an expanding meimjry filter* ^This meanB 
that if the filter is' initialized properly, the estimate^ vill correspond 
to those Qbtalned by designing nonrecursive filters where the data window 

s^^td^at zero and a jpew weight matrix. W(tiT) is computed 4s each new > 

• / . ■■ • ^ .... 

Ts^asurex&ent is made. Obvidpsly as the window expands, the variance of 

.the estimate will- ctecfease; however, if the ^nodel is not exact determinis- 

tic errors will begin to incJrease. 

Iff UsingHhe ei^uatiqpg/'they nmst bcf^ initialized pYoperly if a truly 
• ' ' ' % - 

uhbiased estimate's to be formed. In practice this Is usually done using 
a nonrecuriiive filte?. FoV exact Intt iallzation, driving noise must be 
included in the computation of the nonrecursive filter weights. To mini-- 
mize the computation* the minimum H matrix shojuld be used which oi^ns tlie 
number of rows should ^equal the riumbar of columns. By using this minimum - 
H matrix the filter-heights are computed and "^he initial optimal ^estimate 
\& fonr^d from the actual data* 
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The programs ape written In IX>S FORTRAN. Thoy w^re developed 'on a 
PDF 11/20 with a DOS/BATCH operating system. Printed results are written 
to logical unit 5, which can be assigned at run time to. a line printer, 
CRT terminal » disk file, or other suitable output device. Data Is entered 
from units 6 and- 3,, which can he SRslgned to a/c^rd reader, dink'^data 
fiie^^TTY keyboard^, or any other suitable input device. Plots ciin bo 
pbtained jfith a #T40 graphic^* display terminal and' the plott ins ngibroatines 
provided.. The progtams can be easily modified to use other plotting routinei 



NONRECURSIVE FILTER PROGRAM 

— ■■.■■^■I>>|| (i_aiiiiiiif 1*11 1 ~ " * ' 4 

T^e program for generating the coefficients for non-recursive filters 

is callet^ PROGRAM WINDOW. WINDOW is written in DEC FORTRAN, but may be 

run on other versions of TORTRAN IV with minor n»dif ications. The program 

can call plottixig packages to produce CRT or hardcopy plots of the filter 

^ response to vf^rlotis Inputs . 

Program WINDOW r dads the filter parameters SIGMA, N, M, lA, IPLOT 

• from logical unit 6, where: ' - 

SIGMA determines the input covarlance matrix R. 

If SIGMA > 0.0, R becott^s. 02*1, where a2 
* , - SIGMA and I is the identity matrix. 

If SIGMA < 0.0, the matrix R is read from 
unit 6 in*10F8.2 FORMAT. 

N is the number of points in the window 

(1 < N < 20). ; 

M is the order of the polynomial fit (1 < M < 9). 

I 

lA is the offset a from the first point in the 

window. If lA > 0, the f ijlter predicts tA 

sample times ahead of the most recent sample. 
If lA < 0, the filter .smooths lA sample times 
* behind the oKist recent sample. 

IPLOT is the plotting control variable. If IPLOT 

■ 0, the program finishes after the coeffi- 
cient iMtrlces are computed and printed ^ 
' . out, If IPLOT i 0, an Input signal Is read, 
and the Input points are filtered using the 
^ ^- coefficients computed.' 

• The parameters are re^d in F10.4,.4I3 FORMAT* ' . . . 

' " ... ^ 

^ •^Once the filter parameters have been read, WINDOW generates the re- 

quired S ^Ccovarlance) anfl T matrlopB ,.and computes the weight matrix W- . 

All three matrlcles are then written to logical unit 5. tf IPLOT » 0, 

the program terminates after the weight matrix W has been printed • 

If IPLOT 1* 0, WINDOW reads 101 sample points from' logical unit 3 In 

G15.6 FORMAT.. These points, ^are provided by the user, and are used by 

WINDOW as the filter input signal*. The program filters this Input signal 



using the weight matrix W, and tabulates the input signal, output signal 
(filter response), and error signal (inpirt minus output) for each sample 
time. The tabulated results are printed on unit 5. The sum* of the total 
absolute error is also computed and writteti to unit 5. WtlC)OW then ca^s 
the plotting package routines (subroutine IDIOT) to plot the iupuj&^nd 
output (estimate) signals vs. time. After a PAUSE, the progwfo calls the 
Rlotting package to plot the error signal (input. minus ov»^ut) vs. time. 
^ The plotting packages are included to be used with i^JfTOOW, or WINDOW may 

be changed to call other plotting routines. 

Program WtlJDOW can be modified to write/the derivative estimates. * 

* th / 

Note from equation 2 that the m derivative term is multiplied* by a 

m . 

T /ml fa^or', where T is the sample time. Thus, if the derivative 
estimates are written out, they will be scaled by this factor. To 
illustrate the use of the program, WINDOW was run with the parameters 
SIGMA - 1.0, N - 5r, M - 3, lA - 0, and IPLOT - 0. The filter weighting 
coefficients are listed row by row, and are shown with the input pafa- 
meters and S and T matrices tn Fig. I. The filtef obtained using the 
weight matrix shown estin»tes the input data and the first three deri-. 
vat Ives. Xin equation formsf these estimates are given by 

x[nT] - 0.9857 y{.nT] + 0.05714 y[ (n-l)Tj - 0.0857 y{(n-2)T| 
+ 0.05714 y[(n-3)T] - 0.01429 yI(n-A)T] , , 

c 

X [ nT ] - 1:488 yln Tl - 1.619 yL(n-L)Tj - 7 1 4 ^ LC"r.^J.T 1 



+ 1»048 y {(n -^)Tj - 0.3452 y[(n-4)T ] 

t 



■x[nT] - 0.6429 y[nT] - 1.071 vl(n-l)T] - 0.1429 y{(n~2) Tl 

0.9286 y[(n-3)T] ~ 0. 3571. y [ (n-4)T j « 

T ' 2 
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xInT] - 0.08333 vlnTl - 0.1667"Y[(n-l)TI 

+ 0.1667 y[(n-3)Tf - 0.08333 Y[(n-4 )T1 



RECURSIVE FILTERS 
« 



The Kalman filter prograni Is^ called ADAPT. ADAPT is writtea In DEC 
FORTRAN, but may be ruh on other versions of FORTRAN IV. with minor modi- 
fications. The program can call plotting packages to produce CRT of* 

• • • * 

hardcopy plots of the filter response to various inputs. 

Program ADAPT reads the filter parameters SIQfA, and Q from 
logical unit 6, where 



SIGMA 



M 



/determines the Initial cpvariance 
matrix R. If SIGMA > 0.0^ R becomes 
a2*L, where « SIGMA and I la the 
identity matrix. If SIGMA < 0.x), 
. the matrix R is read from unit 6 In 
10F8.2 FORMAT. 

Is the order of the polynomial fit 
^ (1 < M < 9). , 



Q is the driving noise term. 

The parameters are read in F10.4, I3,F10. A FORMAT. 

Once the filter parameters have been read, ADAPT generates an initial 
S (covariance) and T matrix. ADAPT then generates an initial weight 
matrix W, which is used to inltlaljlze the x vector. The initial covarl- 
ance inatrix is used to initialize the S matrix* 

Once initialized, ApAPT reads 101 sample points from logical unit 3 
in G15.6 FORMAT. These pointy are provided by the user, and are used by 
ADAPT as the Kalman filter input^ The program filters the input using 
equations 85 through 89. The sample number filter \i|put , filter output. 



error signal (input minus output), and KaimaC 
written to unit 5. The total absolute error 



gain are tabulated and 
is also computed and written 
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to unit 5. ADAPT then calls tU^ plotting package routines (subroutine 
IDIO"^ to plot the input and output (estimate) signals* vs. time. 'After 
a PAUSE, the program calls the plotting package to plot the error sig- 
nal (Input minus output) vs. tla». After a second PAUSE, IDIOT Is 



I 



called to plot the Kalman gain vs. time. 

t 

Program ADAPT can be modified to write the derivative estimates. 

th 

Note from equation 2 that the m derivative term is multiplied by a 
m! factor, where T Is the sample tlii». Thus, If . the derivative 

estlmates»|(re written out, they will be scaled by this Tactor. 

/ 

To, Illustrate the use of the program, ADAPT was run with the 
parameters SIGMA - 1.0, M - 3. Q • 0.0. The S (covariance) and T 
matrices used to generate the Initial weight matriK (W) are shown in 
Fig. 2. The S and W matrices are used to Initialize the Kalmmi ftlr 
ters' S and X matrices. The Kalman filter was therv used to filter" 
an ideal sinusoidal signal. A portion of the tabulated results is 
shown in Fig. 3. ' 
GETTING ON LINE 

In order to run progtam WINDOW and ADAPT, first build two files 
named WINDOW. FTN and ADAPT. FTN from the sources provided (card deck 
or paper tape). Also create PLT.FTN and SENDGT.MAC from the sources. 
Compile programs WINDOW and ADAPT with the FORTRAN compiler to create 
WINDOW. OBJ and ADAPT. OBJ. The c?he word Integer option should be 
selected for all compilations. Also, compile PLT.FTN to create 
PLT.OBJ. Assemble SENDGT.MAC under the MACRO aasembler to create 
SENDGT.OBJ. Create a subroutine library called P^il.JLl B.OBJ from 
PLT*OBJ and SENDGT.OBJ' (in that order). Next, LINK WINDOW. OBJ, 
PLTLIB.dfiJ and tl^FTN library to create a file called WINDOW, LDA. 



I- 



/ 



Also., create ADAPT. LDA by LINKlng ADAPT. OBJ, PLTLIB.OBJ and th6 FTN 
library. The files WINDOW. OBJ, ADAPT. OBJ. I'LT.OBJ and SENDGT.OBJ may 
now be deleteci. ^ 
Before rimning WiSDOW or ADAPT, build the file PLOTGT.MAC from 
the source. PLOTGT, is the plotting routine that is loaded into the 
GT40. Assemble and LINK PLOTGT so that PLOTCT.LDA may be loaded Into 
the GT40. After PLOTGT. LDA is running in the GT40, ADAPT. LDA or 
WINDOW. LDA m&y be executed using" a« RUN command. The source listings 
contain additional docun»ntation on these programs. 
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. t**FINITE MEMORV DIGITflL FiLTEff PflCK'ftGE* 

V«RIfl^<CE OF /ERROR* 1 ^©00 I 

SIZE OF WINDOW- 5 I 

CBSJDER OF FIT* 3 , | 

TO PREDICT 0 UNITS ?WftV FROM THE FRONT OF THE WINOOW 

♦ . . . 



M«QIC T MflTRI>< 



1. 


0000 


0. 00@0 


0. 0000 


0 0000 




0000 


-1. e0€» 


1. 0000 I 


-1. 0000 


1 J- 


0000 


~2. 000 


4. 000 


-8- 000 




0000 


-3 000 


9 000 


-27. 00 


i 1 


0000 


' ~4 000 


16. 00 


-64. 00 



VflRIRNCE PtftTRIX, S ^ 



0. 9857 
1 488 
0. 6429 
0. 82:32E~01 



1, 488 
6. 279 
3 869 
0 3972 



^^29 
•869 
^2. 571 
4167 



0. 332:3E-0i 
0. 5972 
0. 4167 
0 6944E-01 



^IVlbSM WEIGHTS. FROM FF^NT TO BfiCK 



0 98157 



1, 4S8. 



i 



i 



ROW i OF W, CORRESPONDING TO THE INPUT SIGftfH. 
0. 5714E-01 -0 Si571E-01 . 0 5714E-01 -0 1429E-01' 

V 

ROW 2 'of W. CORRESPOND iInG TO DERIVRTIVE NUMfeER 1 



-1. €19 



-0, 5714 



1 04^ 



-0 2452 



0 6429 



ROW A 3 C<F W. CORRESPONDING TO DEPIVRTIVE NUME;ER 2 
-1. 071 -0. 1429 e 9286 -0 2571 



ROW 4 OF W.. I20RRESP0NDING TO DERJVfiTI'.'E NUMBER 3 
0. S323E-01 -0.1667- 3576E-06 9 lie? -6 S233E-01 



Ftg. 1 - Program WINDOW Sample Output 



ERIC 



9 



STRUCTURE DIGITftL FILTER PflCKflQE*** 



DR'IVING NOISE= 0 0000 
ORDER OF FIT* 7 
INITIfiL ERROR VRRIRNCE 



X 0000 



MRQIC T MfiTRIN 



1 0000 

1 0000 

1 0000 

j. . 0000 



0. 0000 
~1. 0000 
-2. 000 

-3 00e 



\ 



0 0000 
1. 0000 
4. 000 
9. ?00 



0 0000 

1. 0000 
-8. 00*0 
-2?. 00 



Vf^IR^4CE MflTRIK.. S 
1. 0000 * . 1. 823 i. 0000 

1 S22 . 14' 72 12. 50' . 

0000 12. 50 11. 50 



0» 1667 



2 611 



2 500 



0 1667 
2. 611 
2. 500^ 

0. 5555 



WIW>OW WEIGHTS 



c < 



FROM FRONT TO BflCK 



1. 0090 



ROW 1 OF W. CORRESPOfiDING .TO THE INPUT SIG^WL, 
0. 4128E-05 -0. 2503E-05 0 9527E-0e > 

5. 



. 1 322 

1 0000 

*0 16,67 



ROW 2 OF W.. CORRESPOND imS TO DERIVATIVE NUMBER 1 
-3. 000 1 500 -0 2222 «. 

*'rOW 3 OF CORRESPOND HiG TO DERIVATIVE NUMBER • 2 

-2. S00 2- 000 -0 5000 

t • 

•S- 

ROW 4 OF W. .CORRESPOJIDING TO DER I VRT I ^E , NUMBER 2 
-0 5000 0 5000 • -0. 1667 - , 



Fig. 2 - Program ADAPT Sample Output 



ERIC 



0 



^ TIME 


' OUTPUT 




ESTIMATE 




ERROR 


4. 0100 ^ 


0,' 3894 




0. 2894 




0 


26825-06 


5.000 


0. 4794 




0, 4794 




0 


i627E-04 


6 000 f . 


0 5€4€ 


■ - 


0. 5646 




0 


S460E-04 


7,000 


0 €442 


- 


0 6441 




0 


1017E-02 


•8 


0 7174 




0 7172 




0 


i562E-02 


9. 000 


0 7832 




0 7821 




. 0 


2239E-02 




0. 8415 




. 0 8412 




0 


3122E-02 


11 00 


. 0. 8912 




p. 8908 




0 


4320E-02 


12 00 


^ 0. 9220 




0, 9314 




0 


6002E-02 


12 00 


0 9626 




0 9627 




0 


S295E-02 


14 00 


0. 9855 




0. 9842 




0 


1140E-0:^ 


154 00 , 


0 9975 




0 9959 




0 


1551E~02 


.115 • 00 


0 9996 




^ 9975 


• 


0 


2086E-02 


17 00 


0. 9917 




0 9889 




0. 


2768E-02 


18 00 


0 ^728 




0 9702 




0. 


2620E-02 


19.. 00 


0. 9462 




0, 9416 




0, 


4669E-02 


20 00 


0 9092 




0, 9024 




0 


5928E^02 


21 00 


0. 8622 




0. S55S 




0. 


74S0E-02 


22 00 


0 3085 




0. 7992 




0. 


9229E-02 


2S 00 


•0. 7457 ' 




' ' 0, 7244 




0. 


li30E-01 


24 00 


0, 6755 




0. 6618 




0. 


1267E~01 


25 00 


0. 5985 ' 




0, 5821 




0 


ie25E-0i 


26 00 


0. 5;L55 




0. 4961' 




0 


1927fe-01 


27 00 


0, 4274 




0. 4047 




0 


2272E-01 


28 00 


0 2250 


- 


0 2086 




0 


2640E-01 


29 00 


0 2292 




0 2088 




0. 


2040E-0i 


20 00 


0, 1411 




0. 1064 




0 


2472E-01 


21 00 


0 415eE- 


-01 


0 2265E- 


02 


0 


2922E-01 


22 00 


-0. 5827E- 


■01 


-0 1025 




0^ 


^Jtt|E-01 


-22. 00 


-0 1577 




-0 2070 






BBe>-0i 


24 00 


"0. 2555 




-0. 2100 






MP^-01 


25 00 


-0. 2508 


f 


-0 4106 






i^:lE-0i 


26 00 


-0 4425 




-0. 5077 




6519E-01 


"■7 00 


-0 5298 




-0 6004 




0 


7054E-01 


2i- 00 


-0 6119 




-0. 6876 




0. 


7.576E-0i 


~9' 00 


6378 




-0. 7685 




0. 


S076E-01 




-0 7568 




-0 8422 




0 


8542E-01 


41 00 


-0. S182 




-0. 9080 




0 


89e8E^01 


• 00 . 


-0 8716 




-0 9649 , 




0 


9226E-01 


+ ~ 0«.^ 


-0 9162 




-1 012 




0 


9627E-01 


44 00 


-0, 9516 




-1 050 




0 


9857E~01 


45 00 


-0. 9775 




-1. 077 


i 


■ 0 


9981E-01 


46 00 


-0 9927 




-1 094 




0 


9997E-01 


47 00 


-0. 9999 




~1. 099 




0 


9S91E-01 


4?: 00 


-0. 9962 




-1. 092 




0 • 


9648E-01 


43. 00 . 


-0 9825 




-1. 075 




0 


9255E-01 


5f . 00 


-0. 9589 




-1 046 




0 


8697E-01 


51 ©a 


-0. 9258 




-1. 005 




0. 


79e2E-ei 


52 00 • 


-0. 8825 




* -0 9529 




0 


7029E-01 


52^<»00 

•''54.' 00^^^ 


-0.. 8222 




-0. 8914 




0. 


5915E~01 


-0. 7728 




-0, 8186 




0 


4580E-0i 


55, 00 


■ "0, 7^.55 




-0 7258 




0 


2025E-01 


56, 00 


-0.;e212 




-0, 6427 




0. 


124iE-0i 


57, 00 


5507 




-0 5429 




-0 


7760E--e2 


58. 00 


-0 -4646 




-0 4242 




-0 


2022E-0i 


59, 00 


-0, 2729 




-0, 2186 




-0 


5528E-01 


€0, 00 


"0. 2794 




-0 1968 




-0'3265E-01 


61 00 


-0 1822. 




-0 6975E- 


•01 


-0 


1124 



KflL GI=^IN 
0 9857 

' 0. 9695 
0 94S2 
0 9162 
0 8794 
0 8419 
0 8052 
0. 7707 
0. 7282 
0. 7079 
0, 6797 
0 6524 
0 6289 
0 6061 
0 5848 
0 5649 
0. 5462 
0 5287 
0. 5122 
0. 4968 
0 4822 
0. 4684 
0 45S2 
0 4430 



0 
0 
0 

0 
0 



421: 



4201 
4095 

2995 
2899 
0 2807 
0 2720 
0 2627 
0 2557 
0 2480 
0 2407 
0 2227 
0. 2269 
0 2204 
:i42 
Z.082 
2024 
0 2969 
0 2915 
0 286.2 
0. 2812 
0 2765 
0 2718 
0 2672 
0 2629 
0 25S7 
0 2546 
0 2507 
0 3468 
ii«2421 ' 
0* 2295. 
9 2260; 

2225 
0 2292 



Fig. 3 - Program ADAPT Sample Output 



0.9857 
1.488 
0.6429 
0.08333 



d-. 05714 
v-i.619 
-1.071 
-0.1667 ., 



-0.085? 
-0.5714 
-0.1429 
0.9 • 



\^ "Weight Matrix 



0.05714 
1.048 
.0.9286 
0.1667 




0.9857 


1.488 


0.6429 


0.08;333 


1.488 


6.379 


3.869 


0.5972 


0.6429 


3.869 


2.571 


0.4167 


0.08333 


0.5972 


0.4167 


0.06944 



Covariance Matrix 



Table 1 Optimal Weight Matrix and Covariance Matrix For 

a Five Point Nonrecurslve Filter Using a Third- 
Order Polynomial Model., The Filter Estimates 
the Data at the End of the Window 



-0.08571 
-0.08333 
0.1429 
0.08333 



0.3429 
0.6667 
-0.07143 
-0.1667 



0.4857 
0.0 

•0.1429 
0.0 



0.4857 



.0.0 



-0.1429 



0.0 



Weight Matrix 



0.0 



. 0.9828 



0.0 



-0.2161 



0. 3429 
-0.6667 
-0.07143 

0.1667 



-0.1429 
0.0 

p. 07143 
0.0 



-O.OBS/l 
0.083 33 
0.1429 

-0.08333 



0.0 

-0.2361 
0.0 

0.06944 



Covarlance Matrix 



Table 2 Optimal Weight Matrix and Covariance Matrix For 
, a Five Point Nonrecursive Filter Usfiig a Third 

Order Polynomial Model. The Filter Estlwates 
the Data in the Center of the ii/lndow 
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